package fr.cephb.dbsnptools;


import java.io.BufferedReader;
import java.io.ByteArrayOutputStream;
import java.io.File;
import java.io.FileInputStream;
import java.io.IOException;
import java.io.InputStream;
import java.io.PrintWriter;
import java.io.StringWriter;
import java.util.ArrayList;
import java.util.HashSet;
import java.util.Set;
import java.util.List;
import java.util.concurrent.Callable;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;
import java.util.logging.Logger;
import java.util.regex.Pattern;

import org.xml.sax.InputSource;
import org.xml.sax.SAXException;

import fr.cephb.bio.AcidNucleicUtils;
import fr.cephb.bio.FastaMapping;
import fr.cephb.io.IOUtils;

/**
 * FastMapping of sequences (short reads) using a prefix array and perfect match
 * @author pierre
 *
 */
public class FastMap
	{
	 /** logger */
	private static final int START_FROM_INDEX_0=0;
	/** logger */
    protected static final Logger LOG= Logger.getLogger("fr.cephb");
    /** splitter for SNP context */
    private final Pattern slash=Pattern.compile("[/]");
    /** mapping chromosome /file */
	private FastaMapping mapping=new FastaMapping();
	/** input sources (file,url) containing the reads */
	private Set<String> shortReadFiles=new HashSet<String>();
	/** size of short read */
    private int flanking=60;
    /** max length of the SNP variation */
    private int max_snp_length=100;
    /** maximum number of hits of one SNP on a chromosome */
    private int upper_limit=1000;
    /** sequence of the genomic chromosome */
    private byte chrom_sequence[];
    /** prefix array */
    private int pointers[];
    /** number of threads */
    private int num_threads=2;
    /** number of short read per thread */
    private int num_reads_per_thread=10000;
    
    /**
     * A short read and the information about it (original mapping )
     *
     */
    private static class ShortRead
    	{	
    	String name=null;
    	String context=null;
    	@SuppressWarnings("unused")
		String chromosome=null;
    	Character strand=null;
    	@SuppressWarnings("unused")
    	Integer start=null;
    	@SuppressWarnings("unused")
    	Integer end=null;
    	byte[] seq5=null;
    	String seqMid=null;
    	byte[] seq3=null;
    	}
    
    /** */
    private static abstract class ShortReadsCallable
    	implements Callable<Boolean>
    	{
    	protected String chromName;
    	private List<ShortRead> array;
    	ShortReadsCallable(String chromName,List<ShortRead> array)
    		{
    		this.chromName=chromName;
    		this.array=new ArrayList<ShortRead>(array);
    		}
    	@Override
    	public Boolean call() throws Exception
    		{
    		while(!this.array.isEmpty())
    			{
    			int last=this.array.size()-1;
    			this.map(this.array.get(last));
    			this.array.remove(last);
    			}
    		return true;
    		}
    	public abstract void map(ShortRead shortRead) throws Exception;
    	}
    
    private class MapPairsCallable
	extends ShortReadsCallable
		{
		MapPairsCallable(String chromName,List<ShortRead> array)
			{
			super(chromName,array);
			}
		@Override
		public void map(ShortRead shortRead) throws Exception
			{
			mapPair(this.chromName,shortRead);
			}
		}
    
    
    /**
     * 
     * do something with the current named chromosome and
     * this input source 'f'
     *
     */
    private interface Handler
		{
		public void handle(String chromName,String f) throws IOException;
		}
    
   /** generic run */
    private void run() throws IOException
	    {
	    run(new Handler()
	    	{
			@Override
			public void handle(String chromName, String input) throws IOException
				{
				runFastaInput(chromName, input);
				}
			});
	    }
    
    /** align pairs of short seq */
    private void runPairs() throws IOException
	    {
	    run(new Handler()
	    	{
			@Override
			public void handle(String chromName, String input) throws IOException
				{
				runPairsOfSeq(chromName, input);
				}
			});
	    }
    
    /** align pairs of seq, more info in the input */
    private void runPairsWithContext() throws IOException
	    {
	    run(new Handler()
	    	{
			@Override
			public void handle(String chromName, String input) throws IOException
				{
				runPairsOfSeqWithContext(chromName, input);
				}
			});
	    }
    
    /** align pairs of seq , input from XML-dbsnp */
    private void runDBSNP()  throws IOException
	    {
	    run(new Handler()
	    	{
			@Override
			public void handle(String chromName, String input) throws IOException
				{
				mapDBSNP(chromName, input);
				}
			});
	    }
    
    private void run(Handler handler) throws IOException
    	{
    	//loop over all the chromosomes file
    	for(File chromFile: this.mapping.getFiles())
    		{
    		LOG.info("reading "+chromFile);
    		//get the name for this chromosome
    		String chromName=this.mapping.getName(chromFile);
    		
    		//write genomic sequence in this dynamic byte array
    		ByteArrayOutputStream baos=new ByteArrayOutputStream(150000000);
    		//open the genomiic sequence
    		FileInputStream fin=new FileInputStream(chromFile);
    		
    		//ignore the fasta header
    		int c=fin.read();
    		if(c!='>') throw new IOException("not starting with '>' "+chromFile);
    		///skip fasta header
    		while((c=fin.read())!=-1)
    			{
    			if(c=='\n') break;
    			}
    		//read the fasta sequence
    		byte array[]=new byte[1000000];
    		int nRead;
    		while((nRead=fin.read(array))!=-1)
    			{
    			int j=0;
    			//remove blanks
    			for(int i=0;i< nRead;++i)
    				{
    				if( Character.isWhitespace(array[i])||
    					Character.isDigit(array[i])) continue;
    				array[j]=(byte)Character.toUpperCase(array[i]);
    				if(array[j]=='>') throw new IOException("expected only one sequence in  "+chromFile);
    				j++;
    				}
    			baos.write(array,0,j);
    			}
    		fin.close();
    		//get the genomic sequence as an array of bytes
    		this.chrom_sequence=baos.toByteArray();
    		baos=null;
    		LOG.info("OK in byte array. length= "+chrom_sequence.length);
    		
    		//buildt the array of pointer
    		this.pointers=new int[chrom_sequence.length];
    		for(int i=0;i+this.flanking<= this.pointers.length;++i)
    			{
    			this.pointers[i]=i;
    			}
    		
    		//remove short read containing degenerate bases
    		LOG.info("Removing degenerates");
    		int j=0;
    		for(int i=0;i< this.pointers.length;++i)
    			{
    			boolean ok=true;
    			int k=0;
    			while(k< flanking && this.pointers[i]+k< this.chrom_sequence.length)
    				{
    				byte base=this.chrom_sequence[this.pointers[i]+k];
    				if(!(AcidNucleicUtils.isATGC((char)base)))
    					{
    					ok=false;
    					break;
    					}
    				++k;
    				}
    			if(k==0 || !ok) continue;
    			this.pointers[j]=this.pointers[i];
    			j++;
    			}
    		
    		//compact the array of pointers
    		if(j!=this.pointers.length)
    			{
    			int copy[]=new int[j];
    			System.arraycopy(this.pointers, 0, copy, 0, j);
    			this.pointers=copy;
    			LOG.info("Done. Length is "+this.pointers.length);
    			}
    		
    		LOG.info("OK in int array. Now sorting "+this.pointers.length+" items");
    		//sort the suffix array
    		quicksort();
    		
    		//check sort ok
    		for(int i=1;i< this.pointers.length;++i)
    			{
    			if(compare(i-1, i, FastMap.START_FROM_INDEX_0)>0)
    				{
    				throw new RuntimeException("Not sorted at "+i);
    				}
    			}
    			
    		
    		LOG.info("Done. Sorted");
    		
    		/* echo the first bases
    		for(int i=0;i< 100 && i< this.pointers.length;i++)
    			{
    			int start=pointers[i];
    			StringBuilder x=new StringBuilder(flanking);
    			for(int k= 0;k< flanking && start+k< this.chrom_sequence.length;++k)
    				{
    				x.append((char)this.chrom_sequence[start+k]);
    				}
    			LOG.info("["+i+"] "+x);
    			} */
    		
    		for(String input:this.shortReadFiles)
    			{
    			handler.handle(chromName, input);
    			}
    		}
    	}
    
    private static class Hit
    	{
    	char strand;
    	int position;
		public String toString()
			{
			return "("+strand+") at "+position;
			}
    	}
    
    
    private void runPairsOfSeq(String chromName,String input)
    	throws IOException
	
    	{
    	LOG.info("reading short read file: "+input);
    	String line;
    	BufferedReader in=null;
    	Pattern pat=Pattern.compile("[\t]");
    	//future result, tells where are the available thread
    	Future<?> futures[]=null;
    	//thread manager
    	ExecutorService service=null;
    	//index in futures, where are we looking next
    	int index_in_futures=-1;
    	//the short-reads cuumulation. if size >= num_reads_per_thread, submit to a new thread
    	List<ShortRead> array=null;
    	//if multithread, initialize the variables
    	if(this.num_threads>1)
    		{
    		array=new ArrayList<ShortRead>(this.num_reads_per_thread);
    		futures=new Future[this.num_threads];
    		service=Executors.newFixedThreadPool(futures.length);
    		}
    	try{
	    	 in=IOUtils.openReader(input);
	    	while((line=in.readLine())!=null)
	    		{
	    		if(line.startsWith("#")) continue;
	    		String tokens[]=pat.split(line);
	    		if(tokens.length<3) continue;
	    		ShortRead shortRead=new ShortRead();
	    		shortRead.name=tokens[0];
	    		shortRead.seq5=tokens[1].toUpperCase().getBytes();
	    		shortRead.seq3=tokens[2].toUpperCase().getBytes();
	    		if(service!=null)
	    			{
	    			array.add(shortRead);
	    			
	    			if(array.size()==this.num_reads_per_thread)
	    				{
	    				//find a place for this new thread
	    				while(true)
	    					{
	    					index_in_futures++;
	    					if(index_in_futures>=futures.length)
	    						{
	    						index_in_futures=0;
	    						}
	    					if( futures[index_in_futures]!=null &&
	    						!futures[index_in_futures].isDone())
	    						{
	    						continue;
	    						}
	    					//place was found, submit new thread
	    					futures[index_in_futures]= service.submit(
	    						new MapPairsCallable(chromName, array)
	    						);
	    					break;
	    					}
	    				array.clear();
	    				}
	    			}
	    		/* monothread version */
	    		else
	    			{
	    			mapPair(chromName,shortRead);
	    			}
	    		}
	    	
	    	if(array!=null && !array.isEmpty())
	    		{
	    		new MapPairsCallable(chromName, array).call();
	    		}
    		}
    	catch(IOException err)
    		{
    		throw err;
    		}
    	catch(Exception err)
    		{
    		throw new IOException(err);
    		}
    	finally
    		{
    		if(in!=null) in.close();
    		if(service!=null) service.shutdown();
    		}
    	}
    
    private void runPairsOfSeqWithContext(String chromName,String input)
    throws IOException
	
		{
		LOG.info("reading short read file: "+input);
		String line;
		Pattern pat=Pattern.compile("[\t]");
		BufferedReader in=null;
		Future<?> futures[]=null;
    	ExecutorService service=null;
    	int index_in_futures=-1;
    	List<ShortRead> array=null;
    	if(this.num_threads>1)
    		{
    		array=new ArrayList<ShortRead>(this.num_reads_per_thread);
    		futures=new Future[this.num_threads];
    		service=Executors.newFixedThreadPool(futures.length);
    		}
		
		try
			{
			in=IOUtils.openReader(input);

			while((line=in.readLine())!=null)
				{
				if(line.startsWith("#")) continue;
				String tokens[]=pat.split(line);
				if(tokens.length<9)
					{
					System.err.println("Illegal line: "+line );
					continue;
					}
			
				ShortRead shortRead=new ShortRead();
				shortRead.name=tokens[0];
				shortRead.context=tokens[1];
				shortRead.chromosome=tokens[2];
				shortRead.strand=tokens[3].charAt(0);
				shortRead.start=Integer.parseInt(tokens[4]);
				shortRead.end=Integer.parseInt(tokens[5]);
				shortRead.seq5=tokens[6].toUpperCase().getBytes();
				shortRead.seqMid=tokens[7];
				shortRead.seq3=tokens[8].toUpperCase().getBytes();
				
				if(service!=null)
					{
	    			array.add(shortRead);
	    			
	    			if(array.size()==this.num_reads_per_thread)
	    				{
	    				//find a place for this new thread
	    				while(true)
	    					{
	    					index_in_futures++;
	    					if(index_in_futures>=futures.length)
	    						{
	    						index_in_futures=0;
	    						}
	    					if( futures[index_in_futures]!=null &&
	    						!futures[index_in_futures].isDone())
	    						{
	    						continue;
	    						}
	    					futures[index_in_futures]= service.submit(new MapPairsCallable(chromName, array));
	    					break;
	    					}
	    				array.clear();
	    				}
	    			}
				else //monothreaded
					{
					mapPair(chromName,shortRead);
					}	
				}
			
			if(array!=null && !array.isEmpty())
	    		{
	    		new MapPairsCallable(chromName, array).call();
	    		}
			}
    	catch(IOException err)
    		{
    		err.printStackTrace();
    		throw err;
    		}
    	catch(Exception err)
    		{
    		err.printStackTrace();
    		throw new IOException(err);
    		}
    	finally
    		{
    		if(in!=null) in.close();
    		if(service!=null) service.shutdown();
    		}
		}
    
    
    private void mapDBSNP(final String chromName,String input)
	throws IOException

	{
	InputStream in=IOUtils.openInputStream(input);
	SnpHandler handler=new SnpHandler()
		{
		@Override
		public void found(SnpRecord rec)
			{
			ShortRead shortRead=new ShortRead();
			shortRead.name="rs"+rec.getRsId();
			shortRead.context=rec.getObserved();
			shortRead.seq5=rec.getSeq5(flanking).toUpperCase().getBytes();
			shortRead.seq3=rec.getSeq3(flanking).toUpperCase().getBytes();
			
			mapPair(
					chromName,
					shortRead
					);
			}
		};
	try
		{
		handler.parse(new InputSource(in));
		}
	catch(SAXException err)
		{
		throw new IOException(err);
		}
	in.close();
	}
    
    private void mapPair(
    		String chromName,
    		ShortRead shortRead
    		)
    	{
    	final int max_gap=Math.max(shortRead.seqMid!=null?shortRead.seqMid.length():this.max_snp_length,this.max_snp_length);
    	byte leftSeq[]=shortRead.seq5;
		byte rightSeq[]=shortRead.seq3;
    	List<Hit> leftHits=new ArrayList<Hit>();
		List<Hit> rightHits=new ArrayList<Hit>();
		call(chromName,shortRead.name,leftSeq,leftHits);
		call(chromName,shortRead.name,rightSeq,rightHits);
		if(rightHits.isEmpty() && leftHits.isEmpty()) return;
		
		System.out.println("#c\t"+shortRead.name+"\t"+chromName+"\t"+leftHits.size()+"\t"+rightHits.size());
		if(rightHits.isEmpty() || rightHits.size()>=upper_limit) return;
		if(leftHits.isEmpty() || leftHits.size()>=upper_limit) return;
		
		StringWriter str=new StringWriter();
		PrintWriter w=new PrintWriter(str);
		
		for(Hit left:leftHits)
			{
			for(Hit right:rightHits)
				{

				if(right.strand!=left.strand) continue;
				if(right.strand=='-')
					{
					Hit tmp=right;
					right=left;
					left=tmp;
					}
				
				
				if(left.position>= right.position) continue;
				if(left.position+leftSeq.length+max_gap <  right.position) continue;
				
				
				
				w.print(
					shortRead.name+"\t"+
					chromName+"\t"+left.strand+"\t"+
					(left.position+leftSeq.length)+"\t"+
					right.position+"\t"
					);
				StringBuilder seqGenome=new StringBuilder(max_gap);
				for(int i=left.position+leftSeq.length;i<right.position;++i)
					{
					seqGenome.append((char)chrom_sequence[i]);
					}
				final String seqGenomic=seqGenome.toString();
				w.print(seqGenomic);
				w.println();
				
				
				
				
				if(shortRead.seqMid!=null &&
						shortRead.strand!=null &&
						shortRead.context!=null
						)
					{
					String orientedSeq=seqGenomic;
					
					if(right.strand=='-')
						{
						orientedSeq=AcidNucleicUtils.reverseComplement(orientedSeq);
						}
					if(!orientedSeq.equalsIgnoreCase(shortRead.seqMid))
						{
						w.println("#!d\t"+shortRead.name+"\t"+shortRead.seqMid+"\t"+seqGenomic);
						}
					
					String contexts[]=slash.split(shortRead.context);
					boolean ok=false;
					for(String ctx:contexts)
						{
						if( (ctx.equals("-") && orientedSeq.isEmpty()) ||
							ctx.equalsIgnoreCase(orientedSeq))
							{
							ok=true;
							break;
							}
						}
					if(!ok)
						{
						w.println("#!x\t"+shortRead.name+"\t"+
								shortRead.strand+"\t"+
								shortRead.context+"\t"+
								seqGenomic);
						}
						
					}
				
				}
			}
		w.flush();
		String content=str.toString();
		if(!content.isEmpty())
			{
			System.out.print(str.toString());
			}
    	}
    
    
    private void runFastaInput(String chromName,String input)
    	throws IOException
        {
        int c;
		LOG.info("reading short read file: "+input);
		BufferedReader in=IOUtils.openReader(input);
		ByteArrayOutputStream shortSeq=new ByteArrayOutputStream(flanking);
		StringBuilder shortName=new StringBuilder();
		byte dna[];
		while((c=in.read())!=-1)
			{
			if(c=='>')
				{
				dna=shortSeq.toByteArray();
				if(dna.length!=0)
					{
					call(chromName,shortName.toString(),dna,null);
					}
				
				shortName.setLength(0);
				shortSeq=new ByteArrayOutputStream(flanking);
				while((c=in.read())!=-1)
					{
					if(c=='\n') break;
					shortName.append((char)c);
					}
				}
			else if(Character.isLetter(c))
				{
				shortSeq.write((byte)Character.toUpperCase(c));
				}
			}
		dna=shortSeq.toByteArray();
		if(dna.length!=0)
			{
			shortSeq=new ByteArrayOutputStream(flanking);
			call(chromName,shortName.toString(),dna,null);
			}
		in.close();
		
        }
    
    private void call(String chromName,String shortName,byte shortSeq[],List<Hit> hits)
    	{
	//System.err.println("Calling "+(new String(shortSeq)));
    	if(shortSeq.length==0) return;
    	call(chromName, shortName, shortSeq,'+',hits);
    	String s=new String(shortSeq);
    	s=AcidNucleicUtils.reverseComplement(s);
    	shortSeq=s.getBytes();
    	call(chromName, shortName, shortSeq,'-',hits);
    	}
    
   
    

    private void call(
    	final String chromName,
    	final String shortName,
    	final byte shortSeq[],
    	final char strand,
    	final List<Hit> hits)
    	{
    	int i=lower_bound(0,this.pointers.length,shortSeq);
    	if(i>= this.pointers.length ) return;

	/*
    	for(int k=0;k< flanking && this.pointers[i]+k< chrom_sequence.length;++k)
		{
		System.err.print((char)(chrom_sequence[this.pointers[i]+k]));
		}
	System.err.println();
	*/

	while(i< this.pointers.length &&
			compare(i, shortSeq,START_FROM_INDEX_0)==0)
    		{
    		if(hits==null)
    			{
	    		System.out.println(
	    			shortName+"\t"+
	    			chromName+"\t"+
	    			strand+"\t"+
	    			this.pointers[i]+"\t"+
	    			(this.pointers[i]+shortSeq.length)
	    			);
    			}
    		else if(hits.size()>upper_limit)
    			{
    			return;
    			}
    		else
    			{
    			Hit h=new Hit();
    			h.position=this.pointers[i];
    			h.strand=strand;
    			hits.add(h);
    			}
    		i++;
    		}
    	}
    
    
    /** C+ lower_bound */
	private  int lower_bound(
				int first, int last,
	            byte shortSeq[]
	            )
	    {
		final /* rm later */ int index_in_short_read=START_FROM_INDEX_0;
	    int len = last - first;
	    while (len > 0)
            {
            int half = len / 2;
            int middle = first + half;
            
            
            if (compare(middle,shortSeq,index_in_short_read)<0)
                    {
                    first = middle + 1;
                    len -= half + 1;
                    }
            else
                    {
                    len = half;
                    }
            
            /* we don't have to compare the two short read from the beginning
             * after each while because we know that the strings start
             * with the same prefix
             */
            /*
             * check this LATER
             *
            while(
            	len>0 && 
            	first+len < this.pointers.length &&
            	(this.pointers[first+len]+index_in_short_read) < chrom_sequence.length &&
            	(chrom_sequence[this.pointers[first]+index_in_short_read] ==
            	 chrom_sequence[this.pointers[first+len]+index_in_short_read]) &&
            	index_in_short_read < this.flanking
            	)
            	{
            	index_in_short_read++;
            	}*/
            }
	    return first;
	    }
    
    
    private void quicksort()
    	{
    	quicksort(0,this.pointers.length-1);
    	}
    
    private void quicksort(final int low,final int high)
		{
		if(high<=low) return;
		int i= partition(low,high);
		quicksort(low,i-1);
		quicksort(i+1,high);
		}
	
	 private int partition(int left, int right)
		 {
	     int i = left - 1;
	     int j = right;
	     while (true)
	    	 {
	    	 int index_in_short_read= START_FROM_INDEX_0;
	         while(compare(++i, right,index_in_short_read)<0)      // find item on left to swap
	             {                               // a[right] acts as sentinel
	             //nothing
	             
	             /** TEST THIS LATER
	             while( this.pointers[i]+index_in_short_read < chrom_sequence.length &&
	                 	this.pointers[right]+index_in_short_read < chrom_sequence.length &&
	                 	(chrom_sequence[this.pointers[i]+index_in_short_read] ==
	                 	 chrom_sequence[this.pointers[right]+index_in_short_read]) &&
	                 	index_in_short_read < this.flanking
	                 	)
	                 	{
	                 	index_in_short_read++;
	                 	} */
	             
	             
	             }
	         
	         index_in_short_read= START_FROM_INDEX_0;
	         while (compare(right,--j,index_in_short_read)<0)      // find item on right to swap
	             {
	             if (j == left) break;           // don't go out-of-bounds
	             
	             /** TEST THIS LATER
	             while( this.pointers[right]+index_in_short_read < chrom_sequence.length &&
		                 	this.pointers[j]+index_in_short_read < chrom_sequence.length &&
		                 	(chrom_sequence[this.pointers[right]+index_in_short_read] ==
		                 	 chrom_sequence[this.pointers[j]+index_in_short_read]) &&
		                 	index_in_short_read < this.flanking
		                 	)
		                 	{
		                 	index_in_short_read++;
		                 	}*/
	             
	             
	             }
	         if (i >= j) break;                  // check if pointers cross
	        exchange (i, j);                      // swap two elements into place
	     	}
	     exchange(i, right);                      // swap with partition element
	     return i;
		 }
    
		private int compare(
			final int index1,
			final int index2,
			final int index_in_short_read
			)
			{
			final int base1= this.pointers[index1];
			final int base2= this.pointers[index2];
			
			int i=index_in_short_read,c1,c2;
			while(i<flanking)
				{
				c1=((base1+i)< chrom_sequence.length ?chrom_sequence[base1+i] :-1);
				c2=((base2+i)< chrom_sequence.length ?chrom_sequence[base2+i] :-1);
				if(c1==-1)
					{
					if(c2==-1) return 0;
					return -1;
					}
				else if(c2==-1)
					{
					return 1;
					}
				int diff =c1-c2;
				if(diff!=0) return diff;
				i++;
				}
			return 0;
			}
		
		private int compare(
				final int index1,
				final byte shortSeq[]    ,
				final int index_in_short_read
				)
				{
				final int genome_pos= this.pointers[index1];
				int i=index_in_short_read,c1,c2;
				while(i< shortSeq.length)
					{
					
					c1=((genome_pos+i)< chrom_sequence.length ?chrom_sequence[genome_pos+i] :-1);
					c2=(i< shortSeq.length ?shortSeq[i] :-1);
					if(c1==-1)
						{
						if(c2==-1) return 0;
						return -1;
						}
					else if(c2==-1)
						{
						return 1;
						}
					
					int diff =c1-c2;
					if(diff!=0) return diff;
					i++;
					}
				return 0;
				}
		
		
		
		private void exchange(int i, int j)
			{
			int tmp=this.pointers[i];
			this.pointers[i]=this.pointers[j];
			this.pointers[j]=tmp;
			}
		
		
	public static void main(String[] args) {
		try
			{
			FastMap app=new FastMap();
			int optind=0;
			int mode=0;
			
			while(optind< args.length)
				{
				if(args[optind].equals("-h") ||
				   args[optind].equals("-help") ||
				   args[optind].equals("--help"))
					{
					System.err.println("Options:");
					System.err.println("This program maps the flanking ends of a SNP on a chromosome.");
					System.err.println(" -h help; This screen.");
					System.err.println(" -t <int> number of thread default:"+app.num_threads);
					System.err.println(" -c <file-fasta> add this genomic chromosome fasta file");
					System.err.println(" -f <int> flanking size. default:"+app.flanking);
					System.err.println(" -L <int> max length of a SNP default:"+app.max_snp_length);
					System.err.println(" -H <int> max hit for a flanking end default:"+app.upper_limit);
					System.err.println(" -x <file> read xml mapping file");
					System.err.println(" -m");
					System.err.println("    0: read fasta");
					System.err.println("    1: read name<tab>seq5<tab>seq3");
					System.err.println("    2: read NCBI dbsnp.xml");
					System.err.println("    3: read name<tab>obs<tab>chrom<tab>strand<tab>start<tab>end<tab>seq5<tab>seqSNP<tab>seq3");
					System.err.println("<files for input>");
					return;
					}
				else if(args[optind].equals("-m"))
					{
					mode=Integer.parseInt(args[++optind]);
					}
				else if(args[optind].equals("-t"))
					{
					app.num_threads=Integer.parseInt(args[++optind]);
					}
				else if(args[optind].equals("-x"))
					{
					app.mapping.read(new File(args[++optind]));
					}
				else if(args[optind].equals("-f"))
					{
					app.flanking=Integer.parseInt(args[++optind]);
					}
				else if(args[optind].equals("-L"))
					{
					app.max_snp_length=Integer.parseInt(args[++optind]);
					}
				else if(args[optind].equals("-H"))
					{
					app.upper_limit=Integer.parseInt(args[++optind]);
					}
				else if(args[optind].equals("-c"))
					{
					File f=new File(args[++optind]);
					String name=f.getName();
					if(name.length()>3 && name.endsWith(".fa"))
						{
						name=name.substring(0, name.length()-3);
						}
					if(name.length()>3 && name.startsWith("chr"))
						{
						name=name.substring(3);
						}
					app.mapping.put(name,f);
					}
				else if(args[optind].equals("--"))
					{
					optind++;
					break;
					}
				else if(args[optind].startsWith("-"))
					{
					System.err.println("Unknown option "+args[optind]);
					return;
					}
				else 
					{
					break;
					}
				++optind;
				}
			
            if(optind==args.length)
                    {
                    System.err.println("Illegal number of arguments.");
                    return;
                    }
            else
                    {
                    while(optind< args.length)
                            {
                          	app.shortReadFiles.add(args[optind++]);
                            }
                    }
            if(app.shortReadFiles.isEmpty())
            	{
            	 System.err.println("no short read files");
            	 return;
            	}
            if(app.mapping.getFiles().isEmpty())
            	{
            	System.err.println("no genomic files");
            	return;
            	}
            switch(mode)
            	{
            	case 0:  app.run(); break;
            	case 1: app.runPairs(); break;
            	case 2: app.runDBSNP(); break;
            	case 3: app.runPairsWithContext(); break;
            	default: System.err.println("bad mode "+mode); break;
            	}
			} 
		catch(Throwable err)
			{
			err.printStackTrace();
			}
		}
	}
